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Abstract 


We  have  modified  the  conventional  cell-linked  list  method  to  reduce  the  number  of 
unnecessary  intemuclear  distance  calculations  in  molecular  simulations  of  systems  containing 
many  particles.  In  the  conventional  method,  the  simulation  space  is  partitioned  into  cells  with 
edge  lengths  no  less  than  the  cutoff  distance  of  the  interaction  potential  (rput).  The  atoms  are 
assigned  to  cells  according  to  their  spatial  positions,  and  all  intemuclear  distances  for  atoms 
within  a  cell  and  atoms  in  the  same  and  nearest  neighbor  cells  are  evaluated.  While  this  method 
ensures  that  the  intemuclear  separation  between  all  atom  pairs  within  r<,u,  is  calculated,  it  allows 
for  unnecessary  intemuclear  distance  calculations  between  pairs  that  are  within  the  volume 
encompassing  the  neighbor  cells  but  that  are  separated  by  more  than  r^u,.  The  modified  method 
presented  here  allows  for  reductions  in  the  cell  sizes  and  Ae  number  of  atoms  within  the  volume 
encompassing  the  neighbor  cells.  These  reductions  decrease  the  number  of  atoms  that  are 
outside  of  the  interaction  range  and  the  number  of  unnecessary  intemuclear  distance  calculations 
while  ensuring  that  all  intemuclear  distances  within  the  cutoff  range  are  evaluated.  We  present 
algorithms  to  determine  the  volume  with  the  minimum  number  of  neighbor  cells  as  a  function 
of  cell  size  and  the  identities  of  the  neighboring  cells.  We  also  evaluate  the  serial  performance 
using  the  modified  form  as  functions  of  cell  size  and  particle  density  for  comparison  with  the 
performance  using  the  conventional  cell-linked  list  method. 
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1.  Introduction 


Popular  molecular  simulation  techniques  such  as  molecular  dynamics  or  Monte  Carlo  are 
used  to  study  the  physical  and  chemical  processes  occurring  in  systems  containing  large  numbers 
of  atoms  at  the  atomic  level  (Thompson  1998).  These  methods  require  evaluation  of  either  the 
total  potential  energy  of  a  system  of  N  atoms  (Vxot)  or  the  gradients  of  the  potential  energy.  The 
total  potential  energy  consists  of  terms  that  describe  the  various  interactions  among  the  atoms  in 
the  system.  These  terms  are  usually  functions  of  internal  coordinates,  such  as  intemuclear 
distances  between  two  atoms,  bond  angles  among  three  atoms,  or  torsional  angles  among  four 
atoms.  For  condensed  phase  modeling,  the  total  potential  energy  is  often  described  as  a  sum  of 
two-body  interactions  over  all  atom  pairs.  The  interaction  terms  are  typically  simple  functions  of 
the  intemuclear  distance  ry  between  atoms  i  and  j: 

VT,=si;v(r,).  (1) 

i=l  j>i 


The  evaluation  of  equation  (1)  and  the  gradients  are  usually  the  most  computationally  demanding 
steps  in  a  simulation,  even  if  the  functional  forms  for  V(ry)  are  extremely  simple.  Bmte  force 
evaluation  of  equation  (1)  requires  the  calculation  of  at  least  N  (N- 1)/2  intemuclear  distances.  In 
a  molecular  dynamics  simulation,  each  integration  step  often  requires  the  evaluation  of  equation 
(1)  and  its  gradients  more  than  once  depending  on  the  integration  scheme  that  is  chosen  (Allen 
and  Tildesley  1990).  It  is  clear  that  methods  to  reduce  the  computational  burdens  associated  with 
numerous  evaluations  of  equation  (1)  are  required.  The  most  obvious  recent  approaches  are  to 
modify  the  codes  for  scalable  platforms.  However,  modifications  of  existing  algorithms 
designed  to  reduce  the  computational  burdens  associated  with  evaluation  of  equation  (1)  can  be 
made  to  increase  the  serial  performance  and  exploit  scalable  architectures  to  achieve  enhanced 
performance.  In  this  work,  we  present  a  modification  of  existing  algorithms  that  were  developed 
to  reduce  unnecessary  computations  of  the  intemuclear  distances  for  atom  pairs  used  in  the 
evaluation  of  equation  (1). 
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Common  strategies  to  reduce  the  computational  demands  associated  with  equation  (1) 
include  the  use  of  simple  functions  to  describe  the  pair  interaction  potentials  and  the  assumption 
that  the  interaction  between  two  particles  is  negligible  beyond  a  certain  cutoff  distance,  rcut-  The 
assumption  of  a  cutoff  distance  in  the  interaction  potential  allows  for  a  reduction  in 
computational  time,  since  the  interaction  between  atoms  separated  by  distances  exceeding  rcm 
does  not  need  to  be  calculated.  Unfortunately,  the  easiest  and  most  direct  way  to  determine  the 
set  of  intemuclear  distances  that  are  within  Tout  is  to  evaluate  all  distances  between  all  pairs  and 
eliminate  those  that  exceed  rcut-  This  step  requires  a  potentially  large  number  of  unnecessary 
calculations  and  might  be  the  most  costly  computational  step  in  such  a  simulation. 

The  order  N  method  described  in  the  preceding  paragraph  is  due  to  the  assumption  of  pair 
interaction  potentials  in  equation  (1).  However,  commonly  used  functions  (such  as  Lennard- 
Jones  or  exp-6  forms)  are  too  simple  to  correctly  model  all  of  the  anisotropies  that  exist  in  many 
systems.  Also,  if  chemical  reactions  in  the  condensed  phase  are  being  simulated,  more 
sophisticated  potential  energy  functions  are  required.  Increasingly  complex  potential  energy 
functions  often  use  many  of  the  intemuclear  distances  evaluated  for  equation  (1)  more  than  once 
per  evaluation  of  potential  energy  or  force.  An  example  is  seen  in  the  potential  energy  function 
used  in  the  simulation  of  detonation  (Rice  et  al.  1996).  Li  this  example,  the  function  that 
describes  the  interaction  for  all  atoms  in  the  system  is 

j>i 

where  the  first  set  of  terms  on  the  right-hand  side  of  equation  (2)  (within  the  square  brackets) 
contains  the  intramolecular  interaction  terms  and  includes  many-body  effects.  The  Vvaw  term  in 
equation  (2)  is  a  modified  Lennard-Jones  potential  that  describes  the  intermolecular  interactions. 
The  intra-  and  intermolecular  interaction  terms  have  different  interaction  ranges  and  thus  sample 
different  sets  of  intemuclear  distances  out  of  the  total  set  in  the  system.  The  many-body  term  in 
the  intramolecular  interaction  portion  of  equation  (2)  has  the  form 
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Evaluation  of  this  term  for  a  single  i-j  atom  pair  in  equation  (1)  requires  knowledge  about  the 
remaining  (N-2)  intemuclear  distances.  If  a  bmte  force  calculation  of  the  entire  set  of 
intemuclear  distances  is  performed  for  each  evaluation  of  the  intramolecular  interaction  between 
all  atom  pairs  during  evaluation  of  equation  (1)  using  a  potential  of  the  form  of  equation  (2),  this 
simulation  becomes  order  N^. 


A  reduction  of  uimecessary  calculations  of  intemuclear  distances  can  be  accomplished 
through  the  use  of  the  Verlet  neighbor  list  (Verlet  1967).  This  method  requires  the  constraction 
of  a  list  of  neighbors  for  each  atom.  An  atom’s  neighbors  are  usually  defined  to  be  all  of  the 
atoms  that  are  within  a  distance  slightly  greater  than  the  range  of  the  interaction  potential. 
Information  about  the  neighbors  is  stored  in  arrays.  For  the  duration  of  the  simulation  or  until 
the  lists  are  updated,  each  atom  is  assumed  to  interact  only  with  the  atoms  on  its  neighbor  list. 
The  intemuclear  distances,  interaction  potentials,  and  forces  are  evaluated  for  each  atom  and  its 
neighbors  only.  The  list  may  be  periodically  updated  to  allow  for  the  movement  of  atoms  into  or 
out  of  the  interaction  range.  Bmte  force  constmction  or  update  of  the  list  requires  the  evaluation 
of  all  N(N- 1)/2  intemuclear  distances.  The  method  has  been  shown  to  be  efficient  when  the 
system  contains  a  relatively  small  number  of  atoms  (Allen  and  Tildesley  1990;  Morales,  RuU, 
and  Toxvaerd  1989).  However,  as  the  system  becomes  larger,  the  memory  requirements  for 
maintaining  the  neighbor  lists  become  prohibitive.  Also,  as  the  mobility  of  the  atoms  becomes 
greater,  either  the  frequency  of  lists  updates  must  increase  or  the  cutoff  distance  used  in  the 
definition  of  die  neighbors  must  increase.  Either  of  tiiese  requirements  increases  the 
computational  demands  of  the  Verlet  neighbor  list  method.  The  example  of  the  detonation 
simulation  is  one  such  case  in  which  the  mass  flow  (moving  at  supersonic  speeds)  would  require 
large  neighbor  cutoff  distances  and  frequent  neighbor  list  updates  (Rice  et  al.  1996). 
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Alternative  methods  for  the  efficient  determination  of  the  interacting  neighbors  for  each  atom 
include  grid  or  cell  approaches  (Allen  and  Tildesley  1990;  Boris  1986;  Lambrakos  and  Boris 
1987;  Brage  1993).  These  approaches  partition  the  simulation  space  into  grids  or  cells,  to  which 
the  atoms  are  assigned  by  virtue  of  their  positions  relative  to  the  cells.  Since  each  cell  has  an 
unchanging  set  of  neighboring  cells  that  contain  the  volume  within  the  distance  tcut  of  that  cell, 
an  atom  associated  with  one  of  the  cells  has  as  its  neighbors  those  atoms  assigned  to  the  same  or 
neighboring  cells.  The  implementations  of  these  methods  usually  assign  the  atoms  to  the  cells  at 
each  integration  step.  However,  the  same  considerations  used  for  the  frequency  of  updating  the 
Verlet  neighbor  lists  are  applicable  here.  There  is  some  overhead  associated  with  these  methods, 
and  they  are  preferable  only  for  systems  that  contain  more  than  1,000  atoms  (Allen  and  Tildesley 
1990).  These  methods  substantially  reduce  the  number  of  unnecessary  intemuclear  distance 
calculations  in  evaluating  equation  (1)  but  do  not  completely  eliminate  unnecessary 
computations. 

In  this  work,  we  report  modifications  to  grid-cell  methods  to  further  reduce  the  number  of 
intemuclear  distance  calculations  in  systems  containing  larger  numbers  of  atoms.  The  approach 
we  present  is  a  modification  of  the  conventional  method  of  cell-linked  lists  as  described  in  detail 
by  Allen  and  Tildesley  (1990).  The  results  show  a  dramatic  decrease  in  CPU  requirements  and 
are  amenable  to  parallelization. 

Bruge  and  coworkers  have  already  provided  geometric  and  systolic  parallelization  schemes 
for  conventional  implementations  of  Verlet  neighbor  lists  and  conventional  cell-linked  lists 
(Bragd  1993;  Bragd  and  Fomili  1990a,  1990b).  These  have  shown  significant  decreases  in 
computation  times,  and  we  refer  the  readers  to  such  information.  Our  intent  here  is  to  modify  the 
algorithms  to  accelerate  both  serial  and  scalable  performance.  We  describe  the  modifications 
and  demonstrate  the  performance  on  serial  platforms  in  this  work.  Future  work  will  focus  on 
scalability  and  further  modifications  to  enhance  performance.  We  are  confident  that  some  of  the 
scalable  methods  set  forth  by  Bruge  and  coworkers  (Bruge  1993;  Brugd  and  Fomili  1990a, 
1990b)  will  be  applicable  to  these  algorithms. 
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2.  Method  of  Cell-Linked  Lists 


2.1  Conventional  Method.  The  conventional  method  of  cell-linked  lists  is  well  described 
by  Allen  and  Tildesley  (1990).  We,  like  they,  describe  our  variation  of  the  method  in  two 
dimensions,  but  the  method  can  be  generalized  to  include  three  dimensions.  The  modification 
we  present  is  similar  to  one  suggested  by  Allen  and  Tildesley  (1990) — ^that  the  cell  size  be 
reduced  so  that  no  more  than  one  atom  can  occupy  a  cell. 

In  the  conventional  method  of  cell-linked  lists,  the  simulation  space  is  partitioned  into  cells, 
the  edges  of  which  being  no  smaller  than  the  cutoff  distance  of  the  interaction  potential.  The 
atoms  are  then  assigned  to  the  various  cells  by  virtue  of  their  position  in  the  simulation  space.  A 
linked  list  of  the  atom  indices  is  created  during  the  sorting  procedure.  Also,  at  the  beginning  of  a 
simulation,  an  array  that  contains  a  list  of  cell  neighbors  for  each  cell  is  created.  The  list  remains 
fixed  unless  the  simulation  space  changes  during  the  simulation  (see,  for  example.  Rice  et  al. 
[1996]). 

A  cell  iceii  has  as  its  neighbors  any  cell  that  contains  at  least  one  point  that  is  within  the 
distance  rcut  of  any  point  within  iceu*  Since  the  conventional  method  requires  that  the  edges  of 
each  cell  be  no  smaller  than  rcut,  each  cell  has  eight  nearest  neighbors  (we  are  assuming  periodic 
boundary  conditions  in  both  dimensions  of  our  two-dimensional  example).  These  requirements 
ensure  that  all  atoms  that  are  within  the  interaction  range  of  any  atom  within  iceu  are  assigned  to 
the  eight  nearest-neighbor  cells  of  iceu  or  iceu  itself.  All  atoms  occupying  cells  other  than  these  are 
outside  the  interaction  range  of  any  atom  located  within  iceii-  Figure  1  illustrates  the  division  of  a 
region  of  the  simulation  space  into  cells.  In  this  figure,  both  the  x  and  y  ceU  edges  (denoted  as  lx 
and  ly  hereafter)  equal  rcut-  Evaluation  of  equation  (1)  occurs  through  looping  over  die  cells  using 
the  linked  list  of  atoms  rather  than  accessing  the  atom  indices  sequentially  as  written  in 
equation  (1). 
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This  method  dramatically  reduces  the  number  of  unnecessary  intemuclear  distance 

» 

calculations  that  would  result  from  a  brute  force  calculation  of  all  N(N-l)/2  intemuclear 
distances.  However,  modifications  can  be  made  to  further  reduce  the  number  of  uimecessary 
distance  calculations.  In  the  conventional  method,  the  distances  between  all  atom  pairs  located 
within  the  rectangular  area  of  91xly  are  calculated.  Assuming  the  limiting  case  1^  =  ly  =  rcut,  the 
area  within  which  all  distances  are  calculated  is  9rcut^.  The  area  within  the  cutoff  radius  for  a 
single  atom  is  only  7r*rcut  Thus,  the  traditional  cell-linked  list  method  calculates  distances 
between  all  atom  pairs  within  an  area  that  is  almost  three  times  larger  (or  more,  since  li  >  rcut, 
where  i  =  x  or  y)  than  that  actually  required  for  an  atom.  This  dramatic  difference  is  illustrated 
by  comparing  the  area  within  the  shaded  circle  centered  on  the  atom  labeled ‘T”  with  the  area  for 
the  cell  containing  T  and  its  neighboring  cells  in  Figure  1.  The  shaded  circular  area  illustrates  the 
range  of  interaction  for  atom  T.  Implementation  of  the  conventional  cell-linked  list  for  this 
example  would  result  in  nine  unnecessary  intemuclear  distance  calculations. 

2.2  Modified  Method  of  Cell-Linked  Lists.  The  main  modification  of  the  method  is  in 
the  definition  of  the  sizes  of  the  cells.  By  dividing  the  simulation  space  into  smaller  rectangular 
cells,  each  atom  is  surrounded  by  a  group  of  cells  that  better  approximates  the  area  of  interaction 
for  that  atom.  For  example,  in  Figure  2,  we  have  divided  the  original  rectangular  cells  from 
Figure  1  into  fourths.  The  length  of  each  cell  is  now  l/2rcut.  The  neighboring  cells  to  that 
containing  the  labeled  atom  T  are  the  surrounding  first  and  second  shells  of  cells.  The  area  for 
this  set  of  neighboring  cells  is  6.25rcut^  which  is  approximately  one-third  smaller  than  the 
rectangular  area  that  would  be  considered  in  the  conventional  approach  (see  Figure  1).  Also,  the 
number  of  imnecessary  intemuclear  distance  calculations  has  been  reduced  to  four.  However,  the 
number  of  neighboring  cells  to  that  containing  atom  T  has  increased  from  8  to  24.  Thus,  the  area 
has  been  substantially  reduced,  but  the  number  of  neighboring  cells  has  increased  by  a  factor  of 
3.  There  is  an  increase  in  memory  requirements  associated  with  the  linked  lists  upon  increasing 
the  number  of  neighboring  cells. 
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Figure  1.  Illustration  of  the  Conventional  Cell  Method  in  Two  Dimensions;  Simulation 
Box  Is  Partitioned  Into  3x3  Square  Cells;  Edge  Length  of  Each  Cell  Is  rcut* 
The  Shaded  Circle  Centered  on  Atom  T  Has  Radius  rcut  and  Denotes  the 
Range  of  Interaction  for  Atom  T.  In  This  Method,  the  Eight  Outer  Cells  Are 
Considered  Neighbors  of  the  Central  Cell  That  Contains  Atom  T. 
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Figure  2.  Illustration  of  the  Conventional  Cell  Method  in  Two  Dimensions;  Simulation 
Box  Is  Partitioned  Into  6x6  Square  Cells;  Edge  Length  of  Each  Cell  Is  rcuf 
The  Hatched  Area  Denotes  the  Cells  That  Are  Not  Considered  to  Be  Neighbors 
of  the  Cell  Containing  Atom  T  in  the  Modified  Cell-Linked  List  Method. 
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In  the  simple  examples  shown  in  Figures  1  and  2,  the  division  of  the  original  cell  size  into 
fourths  has  reduced  the  number  of  unnecessary  intemuclear  distance  calculations  from  nine  to 
four,  and  a  further  reduction  in  cell  size  would  probably  not  result  in  additional  savings.  A 
system  whose  atoms  are  arranged  such  that  the  density  is  not  uniform  might  benefit  from  further 
reduction  of  the  cell  sizes  to  the  point  that  the  sphere  of  interaction  of  an  atom  is  closely 
approximated  by  a  set  of  small  rectangular  cells.  Such  an  example  is  given  in  Figure  3,  which 
has  overlaid  the  positions  of  atoms  behind  a  detonation  wave  (Rice  et  al.  1996)  (a  high  dense 
region)  onto  a  grid  of  cells  with  edge  lengths  lx  ”  ly  =  l/2()rcut.  It  is  clear  that  use  of  cells  with  the 
sizes  shown  in  Figures  1  or  2  would  require  many  unnecessary  intemuclear  distance  calculations 
for  the  atomic  arrangement  in  Figure  3. 

Note  that  many  of  the  neighboring  cells  to  that  containing  the  labeled  atom  T  in  Figure  3  are 
empty.  There  is  overhead  associated  with  determining  whether  a  cell  is  occupied.  Also,  more 
memory  is  required  to  maintain  the  cell-linked  list  and  the  neighbor  list,  since  as  the  cell  sizes 
decrease  the  number  of  cells  and  the  number  of  neighboring  cells  for  each  cell  increase.  While 
this  method  reduces  the  unnecessary  distance  calculations,  there  is  a  point  at  which  the  reduction 
in  the  size  of  the  cell  requires  more  computation  in  overhead  than  it  saves  in  eliminating 
unnecessary  distance  calculations.  The  optimum  cell  size  might  vary  from  machine  to  machine 
and  implementation  to  implementation.  Therefore,  it  is  desirable  to  use  an  algorithm  that  allows 
for  the  cell  size  to  be  changed  easily  to  accommodate  portability. 

2.3  Offset  Mapping  Method.  As  cell  sizes  decrease,  memory  requirements  for  storage  of 
neighbor  information  increase  and  are  potentially  a  limitation  on  the  use  of  the  modified  cell- 
linked  list  scheme.  This  problem  can  be  reduced  by  the  determination  of  neighboring  cells 
through  a  list  of  relative  cell  index  offsets,  similar  in  spirit  to  that  presented  in  the  Monotonic 
Logical  Grid  (MLG)  approach  (Boris  1986;  Lambrakos  and  Boris  1987).  After  partitioning  the 
simulation  space  into  cells,  each  cell  is  assigned  a  grid  cell  index  (i,j,k)  that  corresponds  to  its 
location  in  a  Cartesian  reference  frame  (x,y,z).  Figure  4  illustrates  the  two-dimensional  grid 
overlaid  on  the  simulation  box  shown  in  Figure  2.  In  this  example,  the  grid  indices  are  assigned 
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Figure  3.  Simulation  Box  Partitioned  Into  60x60  Square  Cells;  Edge  Length  of  Each  Cell 
Is  l/20rcut.  The  Atoms  That  Are  Illustrated  on  This  Grid  Were  Taken  From 
Results  of  a  Molecular  Dynamics  Simulation  of  Detonation  and  Correspond  to 
the  High  Dense  Region  Behind  the  Detonation  Front  (Rice  et  al.  1996). 
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1  2  3  4  5  6 


Figure  4.  Illustration  of  the  Conventional  Cell  Method  in  Two  Dimensions;  Simulation 
Box  Is  Partitioned  Into  6x6  Square  Cells;  Edge  Length  of  Each  Cell  Is  l/2rcut; 
“Shaped”  Neighbor  Region  (Shaded  Area)  Illustrated. 
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relative  to  cell  (1,1),  located  at  the  lowermost  cell  on  the  left-hand  side  of  the  figure.  The  cell 
that  contains  the  labeled  atom  T  is  located  at  the  fourth  column  (x  direction)  and  the  fourth  row 
(y  direction).  Thus,  the  grid  index  for  this  cell  is  (4,4).  The  set  of  cells  that  are  within  the 
interaction  range  (rcut)  for  all  points  in  cell  (4,4)  consists  of  the  first  and  second  nearest 
neighbors,  each  of  which  has  a  set  of  grid  indices  that  can  be  described  as  relative  offsets  to  (4,4). 
Each  cell  in  the  simulation  box  has  the  same  set  of  relative  grid  index  offsets  as  (4,4).  This  set 
can  be  stored  in  a  relative  offset  array,  which  is  illustrated  for  this  example  in  the  upper  portion 
of  Figure  5. 

Determination  of  the  relative  cell  index  offsets  of  the  neighbors  is  straightforward, 
particularly  if  the  area  encompassing  the  neighbors  is  rectangular.  In  this  example,  the 
rectangular  area  containing  all  neighboring  cells  has  dimensions  of  IXcax  +  li,  where  i  =  x  or  y. 
However,  the  shape  of  the  area  containing  the  neighboring  cells  is  not  limited  to  a  rectangle. 
Further  reductions  in  unnecessary  distance  calculations  can  result  if  the  area  containing  the 
neighboring  cells  resembles  a  circle.  Since  the  set  of  neighbors  must  contain  all  of  the  area 
within  the  interaction  range  of  any  point  within  the  cell,  we  want  the  minimum  set  of  cells  that 
make  up  this  “neighbor  region.”  Rounding  the  comers  of  the  rectangular  neighbor  region  will 
shape  the  neighbor  region  to  approximate  a  circle.  Again,  using  our  simple  example,  we 
illustrate  this  in  the  shaded  portion  of  Figure  4.  The  rounded  comers  represent  the  portions  of 
circles  with  radius  rcm  that  are  centered  on  the  comers  of  the  cell  that  contains  T.  In  this 
example,  the  number  of  cells  containing  the  neighbor  region  is  the  same  as  that  of  the  rectangular 
area.  However,  as  the  cell  sizes  are  reduced,  the  number  of  cells  containing  the  neighbor  region 
will  be  less  than  those  of  the  corresponding  rectangular  area,  and  the  set  of  cells  contained  in  the 
neighbor  range  will  more  closely  approximate  a  circle. 

To  determine  the  minimum  number  of  cells  contained  in  the  neighbor  region,  we  first  assume 
a  rectangular  simulation  box  that  is  larger  than  twice  the  cutoff  radius  in  all  dimensions.  The  box 
is  then  partitioned  into  cells  of  a  desired  size.  At  this  point,  assume  that  the  central  cell  in  this 
box  has  the  grid  index  (0,0).  Only  the  neighbor  cells  contained  in  one  quadrant  of  this  simulation 
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-2,2  -1,2  0,2  1,2  2,2 

-2,1  -1,1  0,1  1,1  2,1 

-2,0  -1,0  0,0  1,0  2,0 


Figure  5.  Geometric  Representation  of  the  Offset  List,  With  the  Relative  Offset  Numbers 
(Upper  Frame),  and  an  Illustration  of  the  Same  Simulation  Box  as  in  Figure  4, 
Surrounded  by  “Ghost”  Cells  (Hatched  Area)  (Lower  Frame).  This  Is  the 
Geometric  Representation  of  the  Mapping  Array.  The  Numbers  Along  the  Left- 
Hand  Side  and  Top  of  the  Figure  Indicate  the  Packing  of  the  Mapping  Array. 
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box  need  to  be  identified,  since  the  remaining  neighboring  cells  that  occupy  the  other  quadrants 
can  be  generated  from  symmetry.  We  consider  the  top  right-hand  quadrant  in  our  illustration. 

The  process  of  identifying  the  neighbor  cells  in  this  quadrant  begins  with  the  calculation  of 
the  range  of  the  cells  along  the  x  axis.  The  grid  index  for  the  furthermost  neighbor  cell  in  the  x 
direction,  xlen,  is  defined  as  xlen  =  floor  (rcut  /  lx+1)*  Each  cell  (ix,0)  that  is  between  cells  (0,0) 
and  (xlen,0)  is  a  neighbor  of  (0,0).  [Due  to  the  symmetry  of  the  system,  all  cells  (jx,0)  that 
include  or  are  between  (0,0)  and  (-xlen,0)  are  also  neighbor  cells.]  We  then  determine  the  range 
in  the  y  direction  as  follows.  For  each  cell  (ix,0)  including  or  between  (0,0)  to  (xlen,0),  the  grid 
index  for  the  furthermost  neighbor  cell  in  the  y  direction  from  cell  (ix,0)  is  defined  as  ylen(ix)  = 
floor  [sqrt(rcut^  -  [(ix- 1)  *  y^)/ly  +  1].  All  cells  that  include  or  are  between  (ix,0)  to  [ix,ylen(ix)] 
are  added  to  the  list  of  neighbor  cells.  The  process  ensures  that  any  cell  whose  lower  left-hand 
comer  is  less  than  rcut  from  the  upper  ri^t  comer  of  the  cell  (0,0)  is  a  neighbor.  A  sample 
FORTRAN  code  for  this  process  is  given  in  the  Appendix. 

In  a  simulation,  when  a  cell  with  grid  index  (i,j)  is  selected,  the  neighboring  cells  are 
identified  by  simply  adding  the  relative  cell  index  offsets  that  are  determined  at  the  begmning  of 
the  simulation  to  the  cell  grid  index  (see  Figure  5  for  the  simple  example  presented  in  this  work). 
If  only  half  of  the  neighbors  are  required  in  the  calculations,  only  the  offsets  in  the  first  and 
second  quadrants  should  be  used,  except  those  from  (-xlen,0)  to  (-1,0).  In  our  example  here,  in 
which  we  have  found  the  neighbors  in  the  upper  right-most  quadrant,  we  may  just  add  the  offsets 
from  (-ix,l)  to  [-ix,  ylen(ix)]  with  ix  =  - 1  to  xlen. 

This  method  as  described  up  to  this  point  is  sufficient  for  determining  neighbors  for  cells  that 
are  far  enough  from  the  edges  of  the  simulation  box  such  that  none  of  the  neighbors  should  be 
minimum  images.  However,  for  cells  on  or  near  the  edge  of  the  simulation  box,  the  method  fails. 
Again,  we  use  our  simple  example  described  in  Figure  2.  hi  Figure  5,  we  have  reproduced  the 
simulation  box  of  Figure  2  and  surrounded  it  with  a  shell  of  “gjiost”  cells  (hatched  area)  that  is 
two  cells  deep  in  both  dimensions.  Overlaying  the  offset  list  (upper  portion  of  Figure  5)  on  cell 
(6,6),  which  is  the  geometric  equivalent  of  simply  adding  the  offset  list  indices  to  (6,6),  would 
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result  in  identifying  as  neighboring  cells  those  with  grid  indices  ranging  from  (4,4)  to  (8,8).  This 
is  clearly  wrong  because  the  cells  with  indices  greater  than  6  are  not  defined.  To  remedy  this,  a 
mapping  array  has  been  developed  to  correctly  map  the  two-dimensional  relative  grid  cell  index 
offsets  to  the  appropriate  set  of  neighbor  cells  while  properly  taking  into  account  the  boundary 
conditions  of  the  simulation.  The  mapping  array  for  the  simple  example  given  in  this  report 
(assuming  periodic  boundary  conditions  in  both  directions)  is  shown  in  the  bottom  portion  of 
Figure  5.  It  is  constructed  using  the  column  and  row  designators  that  border  the  top  and  left-hand 
side  of  the  two-dimensional  array  in  the  bottom  portion  of  Figure  5.  We  illustrate  its  use  as 
follows.  In  this  simple  example,  one  of  the  neighbor  cells  has  relative  cell  index  offset  (0,2). 
Adding  the  relative  cell  index  offset  (0,2)  to  cell  (6,6)  addresses  mapping  array  element  (6,8). 
According  to  the  mapping  scheme,  the  element  (6,8)  of  the  mapping  array  contains  the  grid  ceU 
index  (6,2).  Cell  (6,2),  which  corresponds  to  relative  cell  index  offset  (0,2),  is  the  appropriate 
neighbor  for  cell  (6,6)  according  to  the  periodic  boundary  conditions  established  for  this 
example.  As  for  the  neighbor  list,  this  map  array  increases  in  size  with  decreasing  cell  size. 

By  combining  the  neighbor  offset  list  with  the  mapping  array,  the  computation  of  and 
memory  used  for  storing  the  neighbor  information  are  kept  at  reasonable  levels,  even  for  very 
small  cell  sizes.  There  are  four  major  arrays  associated  with  this  method.  These  are  the  list,  the 
overlay,  the  listhead  (which  contains  the  index  of  the  particle  that  is  used  to  address  the  element 
of  the  linked-list  array)  (Allen  and  Tildesley  1990),  and  the  map  arrays.  The  size  of  the  list  array 
always  equals  the  number  of  atoms.  The  size  of  the  overlay  array  is  proportional  to  the 
interaction  range  divided  by  the  volume  of  the  cell.  The  size  of  the  listhead  array  equals  the 
number  of  cells,  and  the  size  of  the  mapping  array  equals  the  total  number  of  cells  and  ghost 
cells.  When  the  number  of  cells  is  larger  than  the  number  of  atoms,  then  the  listhead  and 
mapping  arrays  require  the  most  memory.  However,  this  method  becomes  inefficient  before  the 
number  of  cells  equals  the  number  of  atoms.  Therefore,  in  any  reasonable  use  of  this  method,  the 
list  array  has  the  largest  memory  requirement 

2.4  Distance  Lists.  As  noted  earlier,  often  more  complex  functions  reuse  information  in 
the  evaluation  of  equation  (1),  such  as  those  systems  that  use  potentials  described  in 
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equations  (2)  and  (3).  This  form  of  interaction  potential  requires  that  the  intemuclear  distances  be 
used  many  times  in  a  single  evaluation  of  the  potential  energy  or  forces.  Recalculations  within  a 
single  integration  step  significantly  increase  the  computer  time  required  for  a  molecular 
dynamics  simulation.  To  overcome  this  problem,  we  implemented  lists  that  contain  information 
about  atom  pairs  for  reuse  in  the  evaluation  of  the  interaction  potential  and  forces.  This  would  be 
unnecessary  for  models  that  assume  pair-additive  interaction  potentials  such  as  the  Lennard- 
Jones  or  exp-6  potentials,  since  the  intemuclear  distances  for  each  pair  are  used  only  one  time  per 
evaluation  of  forces.  But  for  functions  such  as  those  presented  in  equations  (2)  and  (3),  there  are 
several  terms  that  could  benefit  from  storage  of  the  intemuclear  distances.  These  include  the 
exp(-gik)  terms,  the  f(rik)  terms,  corresponding  derivatives  Xjk  and  yac,  distance  r^,  and  the  atom 
index  of  the  neighbor.  This  information  can  be  generated  before  or  during  every  call  to  the 
potential  energy  and  force  subroutine  using  the  linked-list  method  and  neighbor  list.  If  the 
distance  is  within  the  intramolecular  interaction  range,  all  information  that  can  be  reused  is 
calculated  and  stored.  Given  an  atom  pair  ij,  the  stored  information  corresponding  to  that  pair 
can  easily  be  accessed  during  the  evaluation  of  the  potential  energy  and  forces  for  that  pair. 
Since  the  number  of  atom  pairs  can  be  large  compared  to  the  number  of  atoms,  blocking 
techniques  can  be  used  for  the  storage  of  the  atom  pair  information  to  minimize  the  memory 
required  to  store  the  atom  pair  information.  In  a  blocking  method,  the  atom  pair  information  is 
calculated  and  stored  for  only  a  small  number  of  cells  in  the  simulation  space  at  a  time.  The 
potential  for  these  cells  is  calculated,  and  as  the  atom  pair  information  is  no  longer  needed,  it  is 
replaced  by  atom  pair  information  for  other  nearby  cells  that  are  used  next. 

3.  Results 

Although  the  description  of  the  procedure  given  in  the  preceding  section  is  given  in  two- 
dimensional  terms,  the  method  is  tested  herein  in  a  three-dimensional  application.  Six  cubic 
simulation  boxes  that  differ  in  size  have  been  chosen  to  evaluate  this  methodology.  The  six 
simulation  boxes  consist  of  27  (3x3x3),  64  (4x  4x4),  125  (5x5x5),  216  (6x6x6),  343  (7x7x7), 
and  512  (8x8x8)  cubic  cells.  Each  cell  has  edge  lengths  just  greater  than  rent.  The  different 
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simulation  boxes  are  denoted  hereafter  as  simulation  box  3,  4,  5,  6,  7,  and  8,  respectively.  All 
calculations  were  performed  serially  on  an  SGI  Onyx  with  four  195-MHz  RIOOOO  processors 
with  1.5  GB  of  main  memory  and  4  MB  of  secondary  cache  per  processor. 

The  CPU  time  used  to  evaluate  the  intemuclear  distances  using  this  method  as  a  function  of 
system  size  and  particle  density  is  given  in  Table  1.  For  the  evaluation  using  the  cell-linked  list 
methods,  we  report  only  the  times  for  actual  evaluation  and  do  not  include  any  initialization.  The 
initialization,  which  includes  setting  up  the  mapping  array  and  the  relative  cell  offset  list,  is 
relatively  fast  and  is  only  done  once.  The  times  reported  are  the  averages  for  20  separate 
evaluations  of  neighbors,  and  the  timings  include  the  construction  of  the  linked  lists  for  each 
evaluation.  To  check  the  method,  all  atom  pairs  were  calculated  and  compared  to  those 
calculated  through  from  the  brute  force  method.  In  Table  1,  the  variable  Ndiv  denotes  the  number 
of  divisions  along  an  edge  of  the  simulation  box.  For  example,  for  simulation  box  3,  Naiv  =  3 
partitions  each  of  the  three  edges  of  the  box  into  three  sections.  The  simulation  box  has  a  total  of 
27  cells.  This  value  of  Ndiv  corresponds  to  the  conventional  cell-linked  list  method.  Ndiv  =  0 
indicates  that  the  cell-linked  list  method  has  not  been  used,  and  all  N(N-l)/2  intemuclear 
distances  are  calculated.  The  calculations  for  Ndiv  =  0  are  denoted  as  “bmte  force”  calculations. 

It  has  been  established  that  the  conventional  cell-linked  list  metiiod  is  superior  to  the  brute 
force  approach  for  systems  in  which  the  dimensions  are  large  compared  to  the  cutoff  radius  of  the 
potential  (Bmge  1993).  We  have  seen  the  same  result  in  this  study.  Table  1  gives  the  times  for 
evaluation  of  the  intemuclear  distances  as  a  function  of  particle  density  and  Ndiv  for  the  six 
simulation  boxes  and  the  different  methods.  For  simulation  box  3,  the  execution  times  of  the 
conventional  and  modified  cell-linked  list  methods  for  low  densities  are  greater  than  that  of  the 
bmte  force  method.  At  higher  densities,  there  is  a  slight  speed-up  using  the  modified  cell-linked 
list  method  over  the  bmte  force  approach.  Note  that  for  all  densities  for  simulation  box  3,  the 
conventional  method  is  slower  than  the  bmte  force  method.  For  systems  that  are  larger  than 
simulation  box  3,  however,  the  performance  of  the  conventional  and  modified  cell-linked  list 
methods  given  here  are  superior  to  that  of  the  bmte  force  method.  For  the  largest  simulation  box 
(box  8),  there  is  a  90-97%  reduction  in  CPU  time  over  the  bmte  force  method. 
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Table  I.  Time  (in  Milliseconds)  Required  to  Evaluate  Internuclear  Distances  for  Systems  of  Different  Sizes 
and  Particle  Densities 


No.  of  Atoms  Per  Cell  || 

cs 

trj 

Pi 

t/i 

1 

u 

3 

0 

•.n 

X 

to 

X 

fO 

.0 

1 

•i 

§ 

1 

1 

1 

&s 

illii 

so 

ON 

CO 

0 

cs 

wo 

wo 

00 

ON 

wo 

q 

CO 

VO 

s 

ON 

SO 

WO 

ri 

1 

0 

q 

CS 

wo 

w^ 

00 

00 

w^ 

WO 

VO 

II  609 

q 

wo 

q 

I 

0 

q 

q 

wo 

wo 

wo 

On 

WO 

q 

cd 

so 

wo 

ON 

II  0'9£ 

I 

so 

ON 

ON 

cs 

so 

i 

ON 

CO 

cs 

00 

0 

VO^ 

00" 

so 

so 

VO 

SO 

wo 

CO 

so 

00^ 

wo" 

CO 

CO 

00 

wo 

00 

0^ 

wo 

cs 

1 

00" 

On 

0 

i> 

wo 

1 

wo 

VO 

wo 

? 

SO 

00 

q 

so" 

CO 

ON 

CO 

q 

00 

CO 

00 

00 

CO 

5' 

q 

i>" 

On 

s. 

CO 

ON 

ON 

q 

cs" 

p- 

CO 

00 

so" 

00 

p' 

cs 

00" 

p^ 

CO 

§ 

p^ 

cs 

ocT 

0 

00 

9 

cn 

*d 

& 

Tj- 

0 

uo 

wo 

q 

00 

wo 

ON 

On 

WO 

q 

i> 

wo 

d 

wo 

so 

so' 

cs 

1 

0 

q 

cd 

w^ 

0 

On 

WO 

On 

ON 

WO 

vd 

wo 

q 

q 

CO 

1 

0 

00 

cd 

wo 

q 

00 

wo 

ON 

WO 

q 

cd 

wo 

i 

P 

CO 

CO 

.-H 

Tt 

00 

VO 

00 

fHI^ 

tT 

r- 

o' 

r- 

wo 

00 

0 

VO 

p- 

wo 

00 

r-" 

CO 

On 

00 

VO 

<s 

ON 

1 

•M 

0 

ON 

o\ 

wo 

wo 

0 

CS 

CS 

CO 

q 

cs 

1 

VO 

cs 

00" 

T>H 

q 

so 

? 

ON 

00 

CO 

cs 

cs 

SO 

00 

CO 

CO 

00 

0" 

w^ 

ON 

P- 

00 

CO 

cs 

q 

wo" 

CO 

wo" 

CO 

wo 

00 

cs 

ON 

CS 

wo 

5? 

•d 

cs 

0 

q 

CO 

wo 

q 

wo 

to 

00 

CO 

wo 

q 

p- 

On 

CO 

q 

vd 

cs 

1 

0 

00 

§ 

wo 

cs' 

q 

wo 

wo 

00 

cd 

wo 

? 

rp 

vd 

CO 

§ 

Q 

wo 

0 

1> 

CO 

1 

0 

r- 

g 

cs 

cs 

wo 

q 

wo 

cd 

wo 

q 

vd 

CO 

1 

P 

00 

so 

00 

CO 

K 

'd' 

CO 

CO 

<s 

CO 

CO 

On 

§ 

CO 

00 

CO 

4 

X 

rt 

1 

00 

CO 

q 

cs 

s 

K 

d 

q 

00 

ON 

VO 

q 

wo 

CO 

00 

g 

q 

ON 

Tf 

On 

wo 

X 

wo 

0 

1 

1 

.g 

m 

§ 

00 

VO 

q 

S 

ON 

wo 

■irH 

CO 

cs 

g 

00 

q 

sd 

VO 

wo" 

00 

s 

so" 

p- 

wo 

q 

00" 

TP 

q 

cs 

tn 

« 

T3 

so 

0 

ON 

q 

so 

wo 

q 

CO 

CO 

V3 

1 

0 

A 

1 

1 

o3 

q 

wo 

CM 

1 

0 

i 

0 

c4 

wo 

r- 

so 

§ 

q 

cd 

CO 

VO 

ON 

cs 

sd 

CO 

0 

0 

g 

00 

g 

SO 

5? 

00 

§ 

q 

CO 

wo 

wo 

P 

VO 

cs" 

00 

CO 

cs 

so 

ON 

Si 

00 

VO 

q 

CO 

wo 

q 

wo" 

'tt 

0 

q 

CO 

i 

cs" 

s 

q. 

CO 

wo 

co" 

i 

q 

cd' 

fiO 

S5 

00 

0 

j:3 

I 

wo 

so 

00 

so" 

cs 

q 

wd 

vo" 

00 

IS 

so" 

00 

wo 

q 

p^ 

I 

ON 

s 

(9 

-0 

CO 

0 

CO 

q 

wo 

so 

CO 

q 

0 

c4 

cs 

1 

00 

wo 

<s 

1 

0 

P-* 

CO 

00 

wo 

d 

CO 

cs 

d 

cs 

0 

wo 

q 

wo 

cs 

1 

0 

CO 

1 

0 

1> 

CO 

q 

0 

d 

cs 

cs* 

q 

ci 

1 

i 

P 

0 

CO 

5 

so 

CO 

? 

VO 

CO 

CO 

r- 

CO 

0 

wo 

WO 

s 

cs 

00 

r- 

so 

cs 

ON 

CO 

q 

ON 

WO 

ON 

ON 

ON 

r- 

tT 

VO 

0^ 

q 

cs 

q 

1F>-< 

q 

VO 

wo 

q 

cs" 

q 

wo 

CO 

q 

s 

cs" 

cs 

00 

q 

cs 

B 

q 

cs" 

q 

cd 

?5 

<9 

1 

0 

5 

0 

00 

00 

CO 

CO 

es 

q 

wo 

O) 

wo 

cs 

I 

q 

CO 

r- 

1 

cs 

1 

q 

CO 

cs 

0 

VO 

Pi 

cs 

q 

cs 

r-H 

00 

cs 

CO 

1 

wo 

00' 

r- 

1 

1 

vd 

CO 

J 

0 

cs 

pj 

cs 

00 

q 

1 

cs 

sd 

CO 

1 

VO 

cs 

00 

i 

ON 

WO 

1 

1 

P 

so 

so 

cs 

00 

00 

ON 

00 

1 

ON 

CO 

CO 

c3 

CO 

On 

CO 

R 

VO 

CO 

ON 

00 

r- 

so 

so 

00 

tT 

cs 

00 

cs 

wo 

p- 

CO 

p^ 

Cn 

CO 

cs 

wo 

wo 

ON 

P- 

S 

q 

2* 

a 

s 

so 

ON 

cs 

wo 

00 

cs 

a 

a 

00 

cs 

VO 

cs 

00 

cs 

a 

a 

0 

wo 

cs 

wo 

CO 

•8 

I 


■3 


cj 

O 


I 

o 

o 

6 

a 


i 

Q 


ed 


I 

i 

s 


o 

I 


a 

o 

I 


18 


M  X)  U 


Brute  Force  method  (see  text). 

Conventional  cell-linked  list  method  (sec  text). 


Table  1.  Time  (in  Milliseconds)  Required  to  Evaluate  Internuclear  Distances  for  Systems  of  Different  Sizes 
and  Particle  Densities  (continued) 
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®  Brute  Force  method  (see  text). 

®  Conventional  cell-linked  list  method  (see  text). 


Further  comparison  of  the  modified  method  is  made  to  the  conventional  cell-linked  results 
rather  than  those  using  the  brute  force  method.  Table  1  provides  a  percent  reduction  in  time 
using  the  modified  cell-linked  list  method  over  the  conventional  method,  and  Figure  6  provides 
an  illustration.  Each  curve  in  each  frame  of  Figure  6  shows  that  the  percent  time  reduction  first 
increases  with  increasing  Ndiv,  then  decreases  as  Ndiv  bjecomes  larger.  The  subsequent  decrease 
in  performance  with  increasing  values  of  Ndiv  becomes  more  pronounced  for  systems  with  low 
particle  densities.  For  example,  the  curves  for  densities  of  27  and  64  particles  per  cell  show  that 
the  modified  cell-linked  list  method  is  much  slower  than  the  conventional  method  at  large  values 
of  Ndiv.  Conversely,  the  percent  time  reduction  at  large  values  of  Ndiv  for  high  densities  (>343 
particles  per  cell)  is  only  slightly  less  than  the  maximum  value,  indicating  further  time  reduction 
does  not  necessarily  occur  with  increased  partitioning  of  the  simulation  space  (reduction  in  cell 
size).  This  effect  suggests  that  although  the  number  of  unnecessary  intemuclear  distance 
calculations  is  decreasing  with  increasing  Ndiv  (see  Table  2),  the  computational  costs  for  the 
overhead  associated  with  using  a  smaller  cells  are  increasing  and  will  eventually  outweigh  the 
savings  realized  from  the  reduced  number  of  intemuclear  distance  calculations. 

4.  Conclusions 

It  is  clear  that  as  advances  in  scalable  architectures  continue,  more  sophisticated  molecular 
simulations  requiring  more  atoms  and  more  complex  interaction  potentials  will  be  attempted.  It 
is  because  of  this  expectation  that  we  have  modified  the  traditional  cell-linked  list  method  to 
reduce  unnecessary  intemuclear  distance  calculations  for  larger  and  more  complex  systems.  We 
have  shown  a  significant  increase  in  speed  of  the  evaluation  of  information  needed  for  a 
molecular  simulation  through  the  reduction  of  unnecessary  intemuclear  distance  calculations. 
Although  we  have  developed  this  algorithm  for  acceleration  on  serial  machines,  future  efforts 
will  invoke  strategies  for  further  increased  performance  on  scalable  architectures. 
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Figure  6.  Percent  Time  Reduction  of  the  Modified  Cell>Linked  List  Method  Over  the  Conventional  Method  as  a 
Function  of  the  Number  of  Divisions  (Ndiv)  Along  Each  Edge  of  the  Simulation  Box. 


Table  2.  Number  of  Unnecessary  Intemuclear  Distance  Calculations  for  Various  System 
Sizes  and  Particle  Number  Densities 


No.  of  Atoms  per  Cell  1 

Ndiv 

27 

64 

125 

216 

343 

512 

1  Simuilatii3n  Box  coasistins  of  3x3x3  cubic  cells 

26,432 

174,460 

684,208 

2,257,896 

5,646,287 

12,940,352 

K9H 

238,924 

1,317,668 

5,009,417 

14,745,300 

37,232,143 

82,604,224 

3“ 

238,924 

1,317,668 

5,009,417 

14,745,300 

37,232,143 

l[lllllllll[IKi2S^i£i8^ 

6 

110,516 

688,676 

2,479,917 

7,580,688 

18,662,645 

42,348,736 

9 

75,223 

411,940 

1,449,075 

4,994,196 

12,701,118 

27,410,252 

12 

49,109 

316,868 

1,388,325 

3,800,304 

9416,071 

21.179,584 

15 

45,128 

252,972 

968,867 

2,834,820 

7,208,556 

15,228,564 

18 

40,100 

203,688 

764,841 

2,169,240 

5,614,910 

13,097,784 

21 

37,120 

176,216 

683,984 

1,848,988 

4,804,972 

11,715,904 

- 

SimulatiiDU  Box  consisting  of  4x4x4  cubic  cells  I 

70,304 

436,860 

1,696,928 

5.594,028 

13.769,740 

31,587,760 

UEBH 

1421,824 

7,949,700 

30,299,072 

227,164,436 

505,266,768 

4' 

558,688 

3,100,036 

11,799,072 

34,709,844 

87,868,020 

194,888,272 

8 

293,332 

1,609,092 

6,111,572 

17,727,060 

45,043,784 

99,467,856 

12 

170,464 

948,712 

3,858,284 

11496,116 

29,810,540 

63.704,884 

16 

132,992 

825,732 

3,224,804 

8,838,660 

21,984,672 

49,283,664 

20 

115,276 

603,780 

2,315,872 

6465.260 

16.714,388 

35,911,728 

24 

93,200 

463,576 

1,770,764 

5,227,860 

13,790,044 

30,060,984 

28 

82,340 

427,036 

1433,444 

4,626,168 

12,218,292 

27,882,672 

1  Srauiatioii  Box  coBsistmg  of  5x5x5  wbic  cells  1 

iigggi 

145,428 

881,372 

3,404,448 

11,180,808 

27,490,362 

62,783,136 

5,548,197 

31,114,628 

118,658,052 

353,305,692 

891,621,013 

1,985.184,864 

5' 

1,083,072 

6,026,628 

22,954,927 

67,537,692 

171,020,888 

379,552,864 

10 

574,333 

3,114,628 

11,919,888 

34,368,192 

87,661,399 

193,184,864 

15 

324,597 

2,122,036 

7,765,753 

22,393,692 

57,732,571 

125,439,420 

20 

267,238 

1,584,228 

6,192,889 

17,090,928 

42,125,062 

95.162,464 

25 

222,038 

1,121,124 

4,431,177 

12492,412 

32,066,215 

69,197,976 

30 

175,870 

922,576 

3,583,698 

9,951,192 

29,202,191 

63,713,352 

35 

168,950 

857,812 

3,317,824 

8.815,208 

24,289,163 

53,251,408 

1  SimxdadmBox  of  6x6x63  cubic  celk  I 

ESHI 

259,752 

1,555,996 

5,988,268 

19,614.300 

48,193,916 

109.755,152 

KHI 

16,743,444 

93,988,580 

358,498,232 

1.068,753,540 

2.696.284,912 

6,005,484,784 

6' 

1,863,096 

10,381,028 

39,560,732 

116,411,268 

294,833,524 

12 

995,060 

5,349,092 

20476,920 

59,094,372 

151,113,168 

332,557,552 

18 

644,208 

3,646,756 

13.406,952 

38,402,436 

99,192,204 

216,103,172 

24 

476,948 

2,704,100 

10.576,944 

29,333,220 

73,346,640 

163,167,472 

30 

391,924 

1,870,044 

7,549,532 

21,371,448 

54,740,844 

118,474,024 

36 

299,248 

1,687,464 

6,011,796 

18,939,108 

50,250,076 

108,856.920 

42 

282,396 

1,579,524 

5,666,484 

15,456,412 

41,272,636 

90,945,600 

*  The  number  of  intemuclear  distances  that  are  within  the  cutoff  distance  and  are  required  to  be  calculated  in  an  evaluation  of 


equation  (1)  and  its  derivatives. 

^  Corresponds  to  the  brute  force  evaluation  of  the  N*(N- 1)/2  intemuclear  distances  in  a  system  of  N  particles. 
®  Corresponds  to  the  conventional  method  of  celMinked  lists. 
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Table  2.  Number  of  Unnecessary  Internuclear  Distance  Calculations  for  Various  System 
Sizes  and  Particle  Number  Densities  (continued) 


No.  of  Atoms  per  Cell 

27  64  125  216  343  512 

ESSi 

427,310 

2,508,732 

9,629,888 

31,381,752 

77,230,911 

175,732,480 

42,451,120 

238,425,444 

909,481,487 

2,713,097,076 

6,843,353,865 

15,244,669,440 

T 

2,943,694 

62,700,237 

mimmm 

467,483,959 

mmmm 

14 

1,577,899 

8,456,292 

32,650,934 

93,604,704 

239,630,154 

526,643,712 

21 

1,008,145 

5,767,564 

21,274,439 

60,746,676 

156,949,686 

342,409,612" 

28 

753,414 

4,255,620 

16,652,561 

MIMr.Vg!rM 

35 

612,874 

3,008,764 

11,865,187 

33,597,004 

90,128,832 

196,734,024 

42 

460,190 

2,671,232 

10,396,442 

29,836,104 

78,872,346 

171,449,640 

49 

435,242 

2,468,668 

9,364,454 

mmymm 

143,171,928 

II  ;  ,  ,  '  .'T  litiw  ipif4k  j| 

639,300 

3,787,580 

14,510,808 

47,190,180 

116,163,092 

263,943,792 

■SHI 

94,905,276 

533,066,948 

2,033,457,192 

6,068,049,756 

15,304,238,828 

34,095,663,504 

8' 

4,392,636 

24,507,588 

93,457,192 

275,240,796 

1,547,864,464 

16 

2,368,072 

12,580,036 

48,710,196 

139,378,524 

357,232,184 

784,501,136 

24 

1,503,420 

8,582,740 

mihkmm 

510,266,172 

32 

1,134,968 

6,309,060 

24,695332 

69,107,820 

mkkummm 

40 

919,568 

4,477,816 

17,572,392 

49,688,288 

133,800,800 

292,228,740  || 

48 

693,232 

15,616,800 

56 

695,396 

3,785,284 

13,906,096 

35,478,860 

95,880,684 

211,718,836  II 

*  The  number  of  intemuclear  distances  that  are  within  the  cutoff  distance  and  are  required  to  be  calculated  in  an  evaluation  of 


equation  (1)  and  its  derivatives. 

^  Corresponds  to  the  brute  force  evaluation  of  the  N*(N- 1)/2  intemuclear  distances  in  a  system  of  N  particles. 
®  Corresponds  to  the  conventional  method  of  cell-linked  lists. 
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!  Cutoffr  is  the  cut  off  radius.  Maxdim  is  the  maximum  coordinates  for  the  simulations, 

!  mindim  is  the  minimum.  Ndiv  is  the  number  of  divisions  that  the  simulation  is  divided  into. 

!  All  of  these  are  arrays  of  length  2. 
cr2  =  cutofff  *  cutofff 
clen  =  (maxdim  -  mindim)  /  ndiv 
len  =  int(cutofff  /  clen)  +  1 
maxlen  =  len(2) 

!  Iterate  from  the  cell  immediately  next  to  the  test  cell  to  the  last  cell  in  the  x  direction. 

!  Since  the  height  above  the  test  cell  is  always  the  same  as  the  height  above  the  cell 
immediately 

!  next  to  it  we  don’t  calculate  it  here.  We  start  at  2  just  for  array  index  reasons, 
do  i  =  2,len(2)  +  1 

!  Calculate  the  height  above  the  current  cell. 

lengths(i)  =  floor(sqrt(cr2  -  ((i  -  2)  *  clen(2))**2)/clen(l)  +  1) 

enddo 

!  Taking  advantage  of  the  above  mentioned  symmetry 
lengths(l)  =  lengths(2) 
n  =  0 


!  Time  to  replicate  the  cells  for  all  quadrants  and  create  the  offset  list. 
!  Loop  over  every  cell  along  die  x  dimension, 
do  i  =  -len(2),len(2) 

ai  =  abs(i)  +  1 

!  Loop  over  every  cell  along  the  y  axis  for  column  i 
do  j  =  -lengths(ai),lengths(ai) 

!  Don’t  include  cell  (0,0) 
if(i  .ne.  0  .or.  j  .ne.  0)  then 

n  =  n+  1 
overlay(n,l)  =  i 
overlay(n,2)  =  j 

endif 

enddo 

enddo 


! 
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